theta = linspace(0, 2*pi, 30);

% Define the amplification factor for Beam-Warming method
% For a >= 0 case (equation 12.90a)
G_pos = @(mu, theta) 1 - (mu/2)*(3 - 4*exp(-1i*theta) + exp(-2i*theta)) ...
                     + (mu^2/2)*(1 - 2*exp(-1i*theta) + exp(-2i*theta));

% Create figure
figure('Position', [100, 100, 1000, 800]);

% Plot for mu = 0.8 (a >= 0 case)
subplot(2,2,1);
mu = 0.8;
G = G_pos(mu, theta);
plot(real(G), imag(G), 'bo', 'MarkerFaceColor', 'b', 'MarkerSize', 6);
hold on;
plot(cos(theta), sin(theta), 'r--'); % Unit circle
axis equal;
axis([-1.5 1.5 -1.5 1.5]);
title(['\mu = ', num2str(mu)], 'FontSize', 12);
xlabel('Re(G)');
ylabel('Im(G)');
grid on;
legend('Beam-Warming', 'Unit Circle', 'Location', 'southwest');

% Plot for mu = 1.6 (a >= 0 case)
subplot(2,2,2);
mu = 1.6;
G = G_pos(mu, theta);
plot(real(G), imag(G), 'bo', 'MarkerFaceColor', 'b', 'MarkerSize', 6);
hold on;
plot(cos(theta), sin(theta), 'r--'); % Unit circle
axis equal;
axis([-1.5 1.5 -1.5 1.5]);
title(['\mu = ', num2str(mu)], 'FontSize', 12);
xlabel('Re(G)');
ylabel('Im(G)');
grid on;

% Plot for mu = 2.0 (a >= 0 case)
subplot(2,2,3);
mu = 2.0;
G = G_pos(mu, theta);
plot(real(G), imag(G), 'bo', 'MarkerFaceColor', 'b', 'MarkerSize', 6);
hold on;
plot(cos(theta), sin(theta), 'r--'); % Unit circle
axis equal;
axis([-1.5 1.5 -1.5 1.5]);
title(['\mu = ', num2str(mu)], 'FontSize', 12);
xlabel('Re(G)');
ylabel('Im(G)');
grid on;

% Plot for mu = 2.4 (a >= 0 case)
subplot(2,2,4);
mu = 2.4;
G = G_pos(mu, theta);
plot(real(G), imag(G), 'bo', 'MarkerFaceColor', 'b', 'MarkerSize', 6);
hold on;
plot(cos(theta), sin(theta), 'r--'); % Unit circle
axis equal;
axis([-1.5 1.5 -1.5 1.5]);
title(['\mu = ', num2str(mu)], 'FontSize', 12);
xlabel('Re(G)');
ylabel('Im(G)');
grid on;

sgtitle('Absolute Stability of Beam-Warming Method (a \geq 0 case)', 'FontSize', 14);